# generate maps
# -------------------------------------------------------------------- #
# prelims ####
# -------------------------------------------------------------------- #
# clear environment
rm(list = ls())

# import relevant libraries
library(haven)
library(rgdal)
library(raster)
library(readxl) 
library(classInt)
library(RColorBrewer)
library(tmap) 
library(Hmisc) 
library(GISTools) 
library(viridis) 
library(viridisLite) 
library(tidyverse)
library(ggplot2)
library(units)
library(reldist) 
library(sf) 
library(patchwork)

# directories
dir <- paste0("PATH_TO_MAIN_DIRECTORY_HERE/")
dirmaps <- paste0(dir, "figs/maps/")

# -------------------------------------------------------------------- #
# Figure 1: lng barcelona ####
# -------------------------------------------------------------------- #

# lng estimates
filename <- paste0(dir, "data/int/lng_bcn.dta")
LNG <- read_dta(filename) %>%
  subset(Year == 2019) %>%
  select(PlotCode, VLNG100, VLNG200, VLNG350, VLNG500)

# shape
gis_path <- paste0(dir,"orig/catastro/08_900_UH_2020-04-24_SHF/PARCELA/")
PlotShp <- st_read(stringsAsFactors = FALSE,
                   dsn = path.expand(sprintf("%s", gis_path, "PARCELA")), 
                   query="SELECT REFCAT FROM \"PARCELA\" WHERE FECHABAJA >= 30000000 AND PARCELA != '09000'"
) %>%
  st_transform('+proj=utm +zone=31 +ellps=GRS80 +units=m +no_defs')

# merge with shapefile
PlotsMerged <- merge(x = LNG,
                     y = PlotShp,
                     by.x = c("PlotCode"), by.y = c("REFCAT"), all.x = FALSE, all.y = TRUE) %>%
  st_as_sf() 

# barcelona districts
censusfile <- paste0("0301040100_Districtes_ADM_ETRS89.shp")
cen_path <- paste0(dir, "orig/aj_barcelona/0301040100_BCN_UNITATS_ADM/", censusfile)

# shape
BCNDistricts <- st_read(dsn = path.expand(sprintf("%s", cen_path, censusfile)))  %>%
  st_transform('+proj=utm +zone=31 +ellps=GRS80 +units=m +no_defs')

## generate map ####
# -------------------------------------------------------------------- #

scale <- c(0, 0.1, 0.2,  0.3, 0.4, 0.5)

## produce maps
# 100m
LNG100 <- 
  ggplot(PlotsMerged) +
  # lng
  geom_sf(aes(fill = VLNG100), color = NA) +
  scale_fill_viridis("LNG", begin = 0, end = 1, scale, limits=c(0, 0.5)) +
  labs(x="", y="", title="r = 100") + 
  # decrease size of legend
  guides(shape = guide_legend(override.aes = list(size = 0.5))) +
  # decrease size of legend color
  guides(color = guide_legend(override.aes = list(size = 0.5))) +
  # district boundaries
  geom_sf(data = BCNDistricts$geometry, aes(), color = "black", size=10, fill=NA) +
  # theme
  theme(legend.title = element_text(size = 8), 
        legend.text = element_text(size = 8), 
        title =element_text(size=9),
        panel.background = element_rect(fill = "white"), 
        plot.margin = margin(0, 0, 0, 0, "pt"), 
        axis.ticks.y = element_blank(),axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), 
        plot.title = element_text(lineheight=.8, vjust=1), 
        legend.direction="vertical", legend.position = "bottom")

ggsave(filename = paste0(dirmaps,"LNG100_Joint.pdf"), plot = LNG100)

# 200m
LNG200 <- ggplot(PlotsMerged) +
  geom_sf(aes(fill = VLNG200), color = NA) +
  scale_fill_viridis("LNG", begin = 0, end = 1, scale, limits=c(0, 0.5)) +
  labs(x="", y="", title="r = 200") +
  # decrease size of legend
  guides(shape = guide_legend(override.aes = list(size = 0.5))) +
  # decrease size of legend color
  guides(color = guide_legend(override.aes = list(size = 0.5))) +
  # district boundaries
  geom_sf(data = BCNDistricts$geometry, aes(), color = "black", size=10, fill=NA) +
  # theme
  theme(legend.position = "none",
        title =element_text(size=9),
        panel.background = element_rect(fill = "white"), 
        plot.margin = margin(0, 0, 0, 0, "pt"), 
        axis.ticks.y = element_blank(),axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), 
        plot.title = element_text(lineheight=.8, vjust=1))

ggsave(filename = paste0(dirmaps,"LNG200_Joint.pdf"), plot = LNG200)

# 350m
LNG350 <- ggplot(PlotsMerged) +
  geom_sf(aes(fill = VLNG350), color = NA) +
  scale_fill_viridis("LNG", begin = 0, end = 1, scale, limits=c(0, 0.5)) +
  labs(x="", y="", title="r = 350") + 
  # decrease size of legend
  guides(shape = guide_legend(override.aes = list(size = 0.5))) +
  # decrease size of legend color
  guides(color = guide_legend(override.aes = list(size = 0.5))) +
  # district boundaries
  geom_sf(data = BCNDistricts$geometry, aes(), color = "black", size=10, fill=NA) +
  # theme
  theme(legend.position = "none",
        title =element_text(size=9),
        panel.background = element_rect(fill = "white"), 
        plot.margin = margin(0, 0, 0, 0, "pt"), 
        axis.ticks.y = element_blank(),axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), 
        plot.title = element_text(lineheight=.8, vjust=1))

ggsave(filename = paste0(dirmaps,"LNG350_Joint.pdf"), plot = LNG350)

# 500m
LNG500 <- ggplot(PlotsMerged) +
  geom_sf(aes(fill = VLNG500), color = NA) +
  scale_fill_viridis("LNG", begin = 0, end = 1, scale, limits=c(0, 0.5)) +
  labs(x="", y="", title="r = 500") + 
  # decrease size of legend
  guides(shape = guide_legend(override.aes = list(size = 0.5))) +
  # decrease size of legend color
  guides(color = guide_legend(override.aes = list(size = 0.5))) +
  # district boundaries
  geom_sf(data = BCNDistricts$geometry, aes(), color = "black", size=10, fill=NA) +
  # theme
  theme(legend.position = "none",
        title =element_text(size=9),
        panel.background = element_rect(fill = "white"), 
        plot.margin = margin(0, 0, 0, 0, "pt"), 
        axis.ticks.y = element_blank(),axis.text.y = element_blank(), 
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), 
        plot.title = element_text(lineheight=.8, vjust=1))

ggsave(filename = paste0(dirmaps,"LNG500_Joint.pdf"), plot = LNG500)

# combine plots
combined <- LNG100 + LNG200 + LNG350 + LNG500 + plot_layout(guides = "collect") 

# save
ggsave(filename = paste0(dirmaps,"LNG_dist.jpeg"), plot = combined, width = 8, height = 8)

# -------------------------------------------------------------------- #
# Figure 7: new apartment buildings in BCN ####
# -------------------------------------------------------------------- #
## prepare data ####
# -------------------------------------------------------------------- #

# paths
censusfile <- paste0("0301040100_Districtes_ADM_ETRS89.shp")
cen_path <- paste0(dir, "orig/aj_barcelona/0301040100_BCN_UNITATS_ADM/", censusfile)

# districts shape
BCNDistricts <- st_read(dsn = path.expand(sprintf("%s", cen_path, censusfile)))  %>%
  st_transform('+proj=utm +zone=31 +ellps=GRS80 +units=m +no_defs')

# new apartments
load(paste0(dir, "data/int/NewApts_08_900.RData"))

# add buffer
NewApts$Buff350 <- st_buffer(NewApts$cent, dist = 350)

# new apartments 2017-19
NewApts_2019 <- NewApts %>% subset(YearConst == 2019)
NewApts_2018 <- NewApts %>% subset(YearConst == 2018)
NewApts_2017 <- NewApts %>% subset(YearConst == 2017)

## generate map ####
# -------------------------------------------------------------------- #

# name
mapname <- paste0(dir, "figs/maps/NShock_Identification.jpeg")

# districts
map <- 
  ggplot(BCNDistricts) +
  geom_sf(aes()) +
  labs(x= "", y= "") +  
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(),
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), 
        plot.title = element_text(lineheight=.8, face="bold", vjust=1),
        panel.background = element_rect(fill = "white")) 

# add apartments
map_apts <- map +  
  # 2018
  geom_sf(data = NewApts_2019$cent, aes(shape = "2019"), color = "red", show.legend = "point") +
  # 2017
  geom_sf(data = NewApts_2018$cent, aes(shape = "2018"), color = "red", show.legend = "point") +
  # 2016
  geom_sf(data = NewApts_2017$cent, aes(shape = "2017"), color = "red", show.legend = "point") +
  # legend labs
  labs(shape = "Year of Construction")

# add buffers
map_apts_buff <- map_apts + 
  # by year
  geom_sf(data = NewApts_2019$Buff350, color = "orange", fill = NA) +
  geom_sf(data = NewApts_2018$Buff350, color = "orange", fill = NA) +
  geom_sf(data = NewApts_2017$Buff350, color = "orange", fill = NA) +
  # legend position
  theme(legend.position = c(0.9, 0.2))

# store
ggsave(filename = mapname, plot = map_apts_buff, width = 10, height = 10)

# -------------------------------------------------------------------- #
# Figure A3: defining local neighborhoods ####
# -------------------------------------------------------------------- #
## prepare data ####
# -------------------------------------------------------------------- #

# paths
dta_path <- paste0(dir,"data/int/08_900.dta")
gis_path <- paste0(dir, "orig/catastro/08_900_UH_2019-02-06_SHF/PARCELA/")

# numerical data
NumData <- read_dta(dta_path)

# shape
PlotShp <- st_read(dsn = path.expand(sprintf("%s", gis_path, "PARCELA"))) %>%
  subset(PARCELA != 09000) %>% 
  subset(FECHABAJA >= 30000000)

# centroid
PlotShp$cent <- st_centroid(PlotShp)

# generate id variable based on original row allocation
PlotShp$rownum  <- seq.int(nrow(PlotShp))

## polygons for figure ####
# -------------------------------------------------------------------- #

refpol <- 65500 
shapefile <- PlotShp
buffdist <- 100

# choose a polygon
pol <- subset(shapefile, rownum == refpol) 

# generate buffer around it - units are meters
buff <- st_buffer(x = pol, dist = buffdist)

# identify neighboring blocks for each plot
x <- st_intersects(x = buff, y = shapefile$cent) %>%
  unlist() %>%
  as.numeric()

# spatial data containing the polygons
pols_in_buff <- subset(shapefile, shapefile$rownum %in% x) 

# identify reference polygon/s
pols_in_buff$RefPol <- 0
pols_in_buff$RefPol[pols_in_buff$rownum == refpol] <- 1

# keep bigger buffer to plot the map
buffdist2 <- 500
BigBuff <- st_buffer(x = pol, dist = buffdist2)
x <- st_intersects(x = BigBuff, y = shapefile) %>%
  unlist() %>%
  as.numeric()
RestPlot <- subset(shapefile, shapefile$rownum %in% x) 

# smaller buffer
buffdist3 <- 300
BigBuff3 <- st_buffer(x = pol, dist = buffdist3)
x <- st_intersects(x = BigBuff3, y = shapefile) %>%
  unlist() %>%
  as.numeric()
PlotsInMap <- subset(shapefile, shapefile$rownum %in% x) 

# separe refpol and the rest
RefPol <- subset(pols_in_buff, pols_in_buff$RefPol == 1)
RefNeigh <- subset(pols_in_buff, pols_in_buff$RefPol == 0)

## generate map ####
# -------------------------------------------------------------------- #

map <- 
  # this is what will be shown in plot
  tm_shape(PlotsInMap) + tm_borders("black", lwd = .1) + 
  # cartography
  tm_shape(RestPlot) + tm_borders("black", lwd = .1) +
  # color reference polygon
  tm_shape(RefPol) +
  tm_fill(col = "#31a354") +
  tm_borders("black", lwd = .1) +
  # color polygons in buffer
  tm_shape(RefNeigh) +
  tm_fill(col = "#e5f5e0") + # light green
  tm_borders("black", lwd = .1) +
  # add buffer
  tm_shape(buff) +
  tm_borders("red", lwd = .9)

# store
tmap_save(tm = map, filename = paste0(dir, "figs/maps/RelIneqEx.pdf"))

# remove files
rm(map, PlotsInMap, RestPlot, RefPol, shapefile, PlotShp, NumData)

# -------------------------------------------------------------------- #
# closing ####
# -------------------------------------------------------------------- #

rm(list = ls())
